function A=matrixA(term,n)
%term is the # of expansion terms
%n is the # of unknowns in each term

%these are for testing use
% linearpartofbjk=[1,1,5;2,1,4;3,1,3];%zeros(term);
% constpartofbjk =zeros(term);
% constpartofhjk =[1,1,1;1,2,2;2,1,1];%linearpartofbjk;
if(exist('data.mat','file'))
    load data;
else
    disp('Data file reading error!')
end

unknownsposition=(1:n-1)/n;
h=1/n;
coeff1=1/h^2; coeff2=1/2/h;

%each unit matrix has n-1+n-2+n-2 =3n-5 nonzeros
%each unit matrix is n-1 by n-1
%we have term^2 such matrices
nonzeros=term^2*(3*n-5);
ii=zeros(1,nonzeros); jj=zeros(1,nonzeros); kk=zeros(1,nonzeros);
matrixdealt=0;
for row=1:term
    for column=1:term
        base=matrixdealt*(3*n-5);
        bjk=@(x) linearpartofbjk(row,column)*x+constpartofbjk(row,column);
        hjk=@(x) constpartofhjk(row,column)+x-x;
        tmpbjk=bjk(unknownsposition); %length is n-1
        tmphjk=hjk(unknownsposition);
        %diagnonal
        ii(base+1:base+n-1)=(row-1)*(n-1)+1:(row-1)*(n-1)+n-1;
        jj(base+1:base+n-1)=(column-1)*(n-1)+1:(column-1)*(n-1)+n-1;
        kk(base+1:base+n-1)=-2*tmpbjk*coeff1;
        %upper
        ii(base+n:base+2*n-3)=(row-1)*(n-1)+1:(row-1)*(n-1)+n-2;
        jj(base+n:base+2*n-3)=(column-1)*(n-1)+2:(column-1)*(n-1)+n-1;
        kk(base+n:base+2*n-3)=tmpbjk(1:n-2)*coeff1+tmphjk(1:n-2)*coeff2;
        %lower
        ii(base+2*n-2:base+3*n-5)=(row-1)*(n-1)+2:(row-1)*(n-1)+n-1;
        jj(base+2*n-2:base+3*n-5)=(column-1)*(n-1)+1:(column-1)*(n-1)+n-2;
        kk(base+2*n-2:base+3*n-5)=tmpbjk(2:n-1)*coeff1-tmphjk(2:n-1)*coeff2;
        matrixdealt=matrixdealt+1;
    end
end
A=sparse(ii,jj,kk);

    

        





